Fractional transit compartment model for describing drug delayed response to tumors using Mittag-Leffler distribution on age-structured PKPD model

The response of a cell population is often delayed relative to drug injection, and individual cells in a population of cells have a specific age distribution. The application of transit compartment models (TCMs) is a common approach for describing this delay. In this paper, we propose a TCM in which damaged cells caused by a drug are given by a single fractional derivative equation. This model describes the delay as a single equation composed of fractional and ordinary derivatives, instead of a system of ODEs expressed in multiple compartments, applicable to the use of the PK concentration in the model. This model tunes the number of compartments in the existing model and expresses the delay in detail by estimating an appropriate fractional order. We perform model robustness, sensitivity analysis, and change of parameters based on the amount of data. Additionally, we resolve the difficulty in parameter estimation and model simulation using a semigroup property, consisting of a system with a mixture of fractional and ordinary derivatives. This model provides an alternative way to express the delays by estimating an appropriate fractional order without determining the pre-specified number of compartments.


Introduction
Transit compartment models (TCMs) of perturbed tumor growth describe the delay process by which tumors are inhibited by drug administration [1][2][3]. The model generally consists of proliferating and damaged cell compartments. Apoptosis of the damaged cells by drugs does not occur immediately but with a delay. TCMs have been applied to pharmacokinetics and pharmacodynamics (PKPD) to explain the changes in tumors caused by drugs [4][5][6].
Logistic, exponential, and Gompertz models have been widely used to describe cell proliferation. Simeoni et al. propose a new growth model that increases exponentially at the beginning and linearly increases after the threshold [2]. A mortality process was commonly provided with the multiple of first-order degradation and drug concentration given by a pharmacokinetic model, or indirect response [7]. Subsequently, some proliferating cells enter the damaged phases given by the transit processes, and the number of the transit compartments determines the magnitude of delays. Existing TCMs considering transit processes explain cell-phase changes with the same mean residence time (MRT) [8,9]. When the number of the transit compartments is n, the MRT of each phase is considered as 1/k 1 , indicating that k 1 is a degradation rate from one compartment to another (total MRT n/k 1 ). Such common TCMs can be modeled using Erlang distribution, called Erlang TCMs [10]. There are two considerations related to delays in the system of ODEs; the number of compartments (n) of the damaged cells and their degradation rate. Commonly, n is considered as hyperparameter [2,4,11]. This indicates that n is pre-specified for the simulation and then other parameters are estimated using simulation for fitting the experimental data. Our study aims at an alternative way to capture the tumor delay after drug administration. This approach resolves the hyperparameter problem and enables simultaneously estimating parameters to fit the data. For this purpose, we use a system based on a fractional derivative for describing the delay instead of the system of ODEs.
A study for the physical and geometric interpretation of fractional derivatives is shown in [12]. The fractional derivative is used to model biological phenomena and obtain applications in this field. The applications to biological models constructed by fractional-order differential equations produce more realistic results compared to their integer-order counterparts [13,14]. These studies describe that fractional derivatives involve memory and are considerably advantageous for working with biological processes. However, they do not compare integer-order models against experimental data. Therefore, our study constructs a fractional-based TCM model and compared changes in drug delay with ODE-based TCM.
In this paper, we present a fractional TCM based on the Caputo fractional derivative. An age-structured model is first formulated and a Mittag-Leffler (ML) distribution, a non-Markovian distribution, is used for describing the delay [15]. Age is regarded as the time when the cells enter the damaged cells after the drug injection. Consequently, the damaged cells are modeled using an age-structured model and the mortality rate of the damaged cells depends on both age and drug concentration. Mortality rate of damaged cells is formulated using the convolution of the mortality function k out (C, u) of proliferating cells u and a density function f based on age a.
Among density functions, we apply one of the non-Markovian distributions to ML distribution, unlike phase-type distributions including Erlang distribution. The ML distribution is of the form 1 − E α (−t α ), where E a ðtÞ ¼ P 1 n¼0 t n /Γ(1 + α � n) for α 2 (0, 1] [16]. The ML distribution is an exponential distribution for α = 1 and a heavy-tailed distribution for 0 < α < 1. This distribution is related to a fractional derivative that generalizes integer-order derivatives to allow for arbitrary-order derivates [17,18]. We utilize Laplace transform to obtain the simple form because the density function of ML distribution is given by a series. The resulting system is consisting of ordinary and fractional derivative equations that are equations with fractional-order derivatives. In contrast to common TCMs (or Erlang TCMs) given by a system of ODEs to describe the transition of the damaged cells, the fractional transit compartment model (fractional TCM) has a form of a single equation and its mortality is given by a fractional-order derivative. Unlike the pre-specified number of compartments in Erlang TCMs, fractional TCM may resolve the problem of determining the number of compartments and flexibly expresses the delay by estimating an appropriate fractional order. This model is separated into two parts: a delay expressed by the fractional derivative and an ODE equation describing the dynamics of the delayed damaged cells.
The fractional TCM is compared to the Erlang TCM in which the growth rate of proliferating cells is given by Simeoni et al. [2]. To compare them, the simulation of the fractional TCM is carefully considered because fractional TCM has a system of equations with a mixture of ordinary and fractional derivatives in an equation. The property of semigroup of a fractional derivative is applied to the model under some conditions to resolve this challenge [19,20]. Subsequently, model robustness and sensitivity analysis were investigated for model validation. We explored the change of parameter values based on the amount of data, showing that data fitting using fractional TCM is applicable.
The remainder of this paper is organized as follows. In the next section, fractional TCM is derived from the two-compartment model using ML distribution. The PK model, considering dosing regimens, is applied to the model. Model robustness and sensitivity analysis are investigated. This model is compared with Erlang TCM after parameter estimation and particularly, delays are compared with the impacts of the parameter change. In the Discussion section, we summarize the findings and compare the results with previous studies. In the Methods section, we formulate a compartment model of proliferating and damaged cells for describing delayed tumor dynamics.

Derivation of Erlang TCM
The system induced from the method section is summarized as follows: associated with u(0) = u 0 , yðtÞ ¼ R 1 0 �ða; tÞda, and w(t) = u(t) + y(t). Herein, f is a density function depending on age a, and u and y represent proliferating and damaged cells, respectively. k in and k out are the growth and mortality functions of the proliferating cells, respectively. The operator � represents convolution. We assumed that no drug administration indicated that all cells are proliferating, so y(0) = 0. That is, cells begins to be damaged after a drug is administrated. This assumption follows from the Magni et al. [4].
For example, if f is from the point distribution (Dirac-delta function) in Eq (1), that is, f(a) = δ(a − T), then dy dt ¼ k out ðCðtÞ; uðtÞÞ À k out ðCðt À TÞ; uðt À TÞÞ; which represents a delay differential equation (DDE) [21][22][23]. This equation indicates that all individuals had the same residence time T. As the other example, we present an explicit system of ODEs on density function f using Erlang distribution. This system of ODEs is called Erlang TCM that describes the transition process from one compartment to another. Some of tumor cells enter the damaged phases once a drug is administered. If they cascade a multiple-step process with a chain of compartments, it demonstrates the delays motivated by the pathway of signal transduction [8]. We considered an age density of f n in place of f to formulate Erlang TCM from Eq (1), and f n was given by which is a density of Erlang distribution. We have the following relations using the linear chain trick because f n is differentiable [24,25].
If we define E n as E n (t) = (k in � f n )(t)/k 1 , n � 2, then by differentiating E n , we have the following system of ODEs as follows: If y i is the damaged tumor cells with age i, then the total damaged cells y can be considered as y = y 1 + y 2 + � � � + y n . Let y i = E i . Then, Erlang TCM can be derived as . . . dy n dt ¼ k 1 ðy nÀ 1 ðtÞ À y n ðtÞÞ;

Derivation of fractional TCM
A density function f of ML distribution is used for deriving fractional TCM. Since f has the form of series, it is difficult to find the closed form, resulting in analysis intractable and computationally expensive. To resolve these difficulties, Laplace transformation can be applied to Eq (1) to obtain a simple form that results in a form of fractional derivative equation. For this work, we consider a survival function S(t) given by τ > 0 and density function f is given by f(t) = −dS/dt. Taking the Laplace transform of S and f from t to s gives provided with (τs) α < 1. We define a kernel K(t) by To connect K(t) and Laplace transform of (k in � f)(t), we utilize a fractional derivative. The process is summarized as: (i) Define Riemman-Liouville (RL) fractional integral operator of α 2 [0, 1] as for α > 0, and for α = 0, set RL D t a ¼ I, identity operator. Notably, by definition, RL derivative with order 1 satisfies for 0 < α � 1 or for α = 0, set D a t ¼ I, identity operator. By using Taylor expansion for fractional derivative [26], RL and Caputo fractional derivative have the following relationship: Caputo derivative [18] has the form of provided by y(0) = 0. This equation reveals the relation between Laplace transform of Caputo derivative and ordinary derivative. Now let a function G given by and This enables the determination of the Laplace transform of G(t) by yðtÞ as taking the inverse Laplace transform. Thus, fractional TCM is induced as follows.
provided with u(0) = w 0 , y(0) = 0. We restrict the range of α by 0 < α < 1 for the computational purpose because the system returns to the system of ODEs if α = 0 or 1. The second equation in Eq (4) shows that delays of the damaged cells are described by a fractional order and the computation of a fractional-order derivative in time involves an integration over the entire time history of the function (Fig 1(b)). There are many numerical simulation methods for the fractional models, but there are no available methods to implement the system consisting of the mixture of ordinary and fractional derivatives to the best of our knowledge [27][28][29].
To perform the simulation, we need an additional semigroup property and then the Eq (4) is transformed into the simple form to conduct the simulation. We will discuss this property detail in the other subsection.

Growth and mortality functions of the proliferating cells
Simeoni et al. develop the rate change of the proliferating cells u given by a perturbed growth (tumor) function k in (u, w) and mortality function k out (C, u) based on first-order, commonly used in Erlang TCMs [2]. k in and k out are given by For a sufficiently large ϕ (more than 10), k in has approximately exponential growth with a linear rate λ 0 u less than or equal to the threshold w th = λ 1 /λ 0 and nonlinear growth with a rate of λ 1 � u/w otherwise. They provide two growth rates in a single form for computational reasons.

Pharmacokinetic model
The PK model is from the Magni study [4], which also used the growth function given by Simeoni et al., and consists of two compartments, as follows: where q 1 and q 2 are the quantity of drug in the plasma and peripheral compartments, respectively, and V is the volume in plasma. C(t)(ng � ml −1 ) is the concentration, and v(t) is the bolus administration (ng � kg −1 ). The obtained drug concentration C is applied to Erlang and fractional TCMs shown in Eqs (2) and (4).

Parameter values and estimations
The tumor data consists of ten points of tumor size versus time profile, as shown in mouse 150 [4]. The tumor is implanted on day 0 with an initial tumor size w 0 given by 0.0121g, and the drug is administered on day 13 and injected every day for ten days. Their study and PK profile are reproduced in detail, as shown in  fractional TCM [27]. A fractional derivative equation solver can be used after a modification with the semigroup property.

System modification for implementation: Semigroup property
Eq (4) is a fractional system consisting of two equations. One of the equations has ordinary and fractional derivatives. In this case, it is challenging to simulate the system. Particularly, fde12 cannot be directly applied because the system consists of a mixture of ordinary and fractional derivatives in an equation. To resolve this problem, we considered a latent variable z such that z is defined as z ¼ D 1À a t y. Generally, the Caputo derivative dose not satisfy a semigroup property with respect to α of differentiation; however, it possesses a semigroup property under some additional assumptions [19,20], which are D 1 t yðtÞ and D 1À a t ðD a t yðtÞÞ are continuous on R. Then semigroup property is satisfied. The system given by Eq (4) is transformed into a multi-order system of fractionalorder equations and computer simulation such as fde12 can be used. The resulting system is as follows.
where C is from the PK model (Eq (5)). Notably, the damaged cells (y) given by z, are delayed.

Comparison between Erlang and fractional TCMs
Erlang and fractional TCMs given by Eqs (2) and (6)  to determine the optimal number of compartments. Here, N is the amount of data, D n is the nth data, and w n = u n + y n . Furthermore, u n and y n are model outputs of proliferating and total damaged cells at the nth data time. The number is four, which represents four transit compartments, except for the proliferating compartment. Fractional TCM has RMSE = 0.3392 which is comparable to Erlang TCM, which has RMSE = 0.3139 for full data N = 10. From the obtained RMSE, we measure Akaike information criterion (AIC) which is an estimator of prediction error and thereby relative quality of statistical models for data. AIC is given by and k as number of parameters. AIC of Erlang and fractional TCMs are 11.2023 and 14.7538 which are comparable. We note that although the AIC of Erlang TCM is smaller than fractional TCM, the number of compartments for describing delays of the damaged cells is prespecified to simulate the model. This shows that Erlang TCM is necessary for two parameters k 1 and n for capturing the delays, but n is pre-specified through several tests and then k 1 can be estimated. In the case of fractional TCM, however, there are two estimated parameters τ and α to describe the delays, but the estimation simultaneously enables. Moreover, the number of equations is smaller in the fractional TCM. Our approach enables reduce the difficulty to specify n and provides a single equation for the delayed dynamics of the damaged cells.

Model robustness
In the two TCMs, we perform data fitting from partial data for model robustness. Model robustness is that if the outputs or forecasts are consistently accurate even if one or more of the input variables or assumptions are drastically changed due to unforeseen circumstances. We conduct the model robustness test only for the parameters related to the delays. A change of delay is explored by varying two parameters α and τ. We do not investigate the perturbation of η because this is already investigated in Erlang TCM. α and τ are randomly selected from uniform and lognormal distributions. Ten (α, τ) or twenty (both) samples are plotted in Fig 4. The other parameters are the same as we mentioned before. In Fig 4(a), delays are smaller when α is larger. Conversely, delays are larger when τ is larger, as shown in Fig 4(b). In Fig 4  (c), α and τ are simultaneously changed. The results demonstrate that fractional TCM can capture various delays for multiple dosing. Furthermore, the fractional model describes dynamics that tumors begin to increase when dosing is stopped on day 23.

Change of parameter estimation based on the amount of data
We test to find minimal data points to capture the full data set using the model. This work shows how the estimated parameters with fewer data points are close to them obtained by the full data set and proposes minimum data points for the given phenomenon. To that end, we randomly selected two, three, four, five, six and eight points (circle) among ten data points (star), as shown in Fig 5. To compare predictions of the two models, α, η, and τ are estimated   data. Other cases, except for two points, fairly predict full data with less data set. However, Erlang TCM is suppressed over data during the treatment, as shown in (d). We also calculate the change of RMSE and values of the parameters based on the amount of data, as shown in Fig 6. RMSE is smaller when the number of data increases. Moreover, the values of the parameters become stable and the values converged. Thus, parameters estimated by randomly selected six data points are similar to those composed of the full data set, meaning that the fractional TCM fairly captures the full data set.

Sensitivity analysis
Sensitivity analysis is performed to understand parameter influences of the tumor size in fractional TCM. Partial rank correlation coefficient (PRCC) is investigated for global sensitivity. PRCC is the correlation between two variables after removing the effect of variables. Detailed explanation and formula are seen in [30]. PRCC in this study is measured for the relation among tumor sizes and parameters. We select 3000 samples of α, τ, and η using the Latin hypercube sampling (LHS) method to perform PRCC, as shown in Fig 7(a). The obtained sample points of parameters are plotted in Fig 7(b). The values of PRCCs are calculated and plotted over time in Fig 7(c). Before day 13 (no injection), the parameters do not influence tumor size. η has a fully negative influence on the increase of tumor size. τ has a strong positive influence during injections; however, the relation decreases and has negative influence as time elapses. α has a negative influence, but it changes positively. From this analysis, α and τ are important parameters with respect to the change of tumor sizes with drug administration, and we may control these parameters in this model.

Discussion
Responses of cell populations have attracted attention for use in PKPD. TCMs describe delays or aging processes in cell populations after drug administration. TCMs explain delays in which damaged cells are eliminated as transit processes and are usually expressed as an ODE system. We present fractional TCM to describe the delay process. A stochastic process is applied to the transit process by a convolution of the density function and degradation rate. Using Mittag- Leffler distribution and Laplace transformation, we formulate the fractional TCM. This model consists of two equations: proliferating and damaged cells. The latter consists of both ordinary and fractional derivatives that are related to the delay of tumor size after drug administration. Model simulation is performed and compared to Erlang TCM. We utilize the semigroup property for model implementation. Thus, the system consists of three compartments: proliferating cells, delay process given by fractional derivative, and damaged cells. Parameter estimation, robustness, parameter change by the amount of data, and sensitivity analysis are performed to verify the model's quality and usefulness. Similar to Erlang TCM, fractional TCM can be used to describe tumor delay caused by drugs without the pre-specified number of damaged compartments with fewer compartments. This approach has a benefit that enables data estimation simultaneously and alleviates the difficulty of choosing the number of compartments n for the damaged cells.
Our simulation study has some limitations associated with the determination of the time step and fractional order. In simulation implementation, the time step should be sufficiently small. In this study, the time step is 0.0039 (that is, 2 −8 ). Otherwise, model predictions hardly demonstrate data set greater than the given time step, even if parameters are correctly selected. The reason could be the numerical method or fractional order. The calibration of the time step for the simulation is not the main issue in the Erlang TCM but is necessary for the fractional TCM. If the time step is too small, then computational work is expensive, and in the opposite case, the model simulation could be unreliable. Another limitation is fractional order α. In contrast to Erlang TCM which has the parameter k 1 related to mean duration and the number of compartments n, the physical meaning of α is still ambiguous. We utilize α for describing the delay in this study. This study does not determine which models are better; however, fractional TCM may be applicable when the number of compartments for the delay is unknown or data fitting seems unreliable. We present an alternative model for TCM and provide a simple model modification for the implementation of simulation using the semigroup property in the system.
Erlang TCMs are widely used in PKPD [2,4,8], epidemics [31][32][33], and other research areas [11,34,35]. However, the time delay is discretely described with regards to the number of compartments [36]. To describe delays, fractional models are substituted with a fractional derivative instead of the ordinary derivative on the left-hand side of the equation, which explains the abnormal kinetic (power-law kinetic) [37][38][39]. Additionally, the use of the fractional models on the left-hand side is already used in the PK study, demonstrating the fat-tail end behavior of the drugs [40][41][42]. In the SIR model, a model using a stochastic process is proposed to reflect the delay of infected individuals [17,43]. Our model approach proposes a method similar to that of the SIR model to represent the delay of cell apoptosis due to drugs, considering data fitting. Some studies have a single equation with the mixture of ordinary and fractional derivatives, or a system of fractional derivatives only in the equations. However, our system is different because the model is a system consisting of a mixture of fractional and ordinary derivatives in an equation.

Mathematical model formulation: Age-structured-based perturbed tumor model using survival function
We consider a cohort of cells: a group of cells of age in an interval of length Δa. We obtain the Mckendrick-von Foerster equation as follows.
where u = u(t) is the number of proliferating cells with initial u(0) = u 0 and C = C(t) is the drug concentration. The total number of damaged cells at time t is given by yðtÞ ¼ R 1 0 �ða; tÞda. w = w(t) is the total number of cells given by w(t) = u(t)+ y(t). We expected that the unit of age is the same as that of time and so, we assumed da/dt = 1. Additionally, if the mortality rate μ(a, C) depends only on age a, then μ(a, C) = μ(a) [44]. This assumption is reliable because the age depends on the duration after the drug C is injected. At a = 0, the boundary condition is ϕ(0, t) = k out (C, u), indicating that damaged tumor cells of age zero are created by tumor size and drug concentration. The initial condition of ϕ is given by ϕ (a, 0) = 0, indicating that all tumors without drug administration are proliferating cells. Under these conditions, time t is the time since administration of the drug.
The existence and uniqueness of the solution for Eq (7) are proven using the method of characteristics, which can be solved as follows [45].
The mortality rate μ(a) is related to the probability of survival. We consider a stochastic process to verify it. Let S(a) = Pr{T � a} be the probability of survival from zero to age a. S(a)y denotes the cells remaining in the cohort until age a; Sða þ DaÞy À SðaÞy ¼ À mðaÞSðaÞyDa denotes the number of cells that die in Δa. Subsequently, both sides are divided by yΔa and Δa ! 0, dS da ¼ À mðaÞSðaÞ: Because the initial value of S is the probability of survival until age 0, we assume S(0) = 1. By solving this equation, we obtain SðaÞ ¼ expðÀ R a 0 mðaÞdaÞ, whose derivative defines the probability density function f of Pr{T < a}, such that f(a) = −dS/da = μ(a)S(a). This indicates the probability that a cohort dies within age a. Thus, μ can be interpreted as the death-hazard rate. By substituting Eq (8) into Eq (9) with f, we derive the following equation.
associated with u(0) = u 0 , yðtÞ ¼ R 1 0 �ða; tÞda, and w(t) = u(t) + y(t). The mortality (elimination) rate of collection of damaged cells y is delayed and given in the form of convolution.